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Abstract 

A simple algorithm with time complexity 0(M(n) log(n) 2 ) and space complexity 
0(n) for the evaluation of the hypergeometric series with rational coefficients is con- 
structed (M(n) being the complexity of integer multiplication). It is shown that this 
algorithm is suitable in practical informatics for constructive analogues of often used 
constants of analysis. 

Introduction. In this paper we construct an algorithm for the calculation of the 
approximate values of the hypergeometric series with rational coefficients whose imple- 
mentation is simple and which is quasi-linear in time and linear in space on the machine 
Schonhage [1] . Such series are used for the calculation of some mathematical constants 
of analysis and of the values of elementary functions at rational points. 

Sch(FQLIN - TIME I /LIN - SPACE) will be used to denote the class of algo- 
rithms which are computable on Schonhage and are quasi-linear in time and linear in 
space. The main feature of Schonhage is its ability to execute recursive calls of proce- 
dures. Quasi- linear means that the complexity function is bounded by 0(n log (n) k ) for 
some k. 

The main advantage of algorithms based on series expansions is the relative simplic- 
ity of both the algorithms and the analysis of their computational complexity. Besides 
we can compute all the most commonly used constants of analysis using series expan- 
sions. For calculations with a small number of digits after the binary (or decimal) point, 
series are more efficient than other methods because of the small constants in estima- 
tions of their computational complexity. Therefore such algorithms are important in 
computer science for practical applications. 

It is known [2] that linearly convergent hypergeometric series with rational coef- 
ficients can be calculated using the binary splitting method with time complexity 
0(M(n)(log(n)) 2 ) and space complexity 0(nlog(n)) (where M(n) denotes the com- 
plexity of multiplication of re-bit integers). In recent publications, for example [3], 
algorithms based on a modified binary splitting method for the evaluation of linearly 
convergent hypergeometric series with time complexity 0(M(n) (log(n)) 2 ) and space 
complexity 0(n) are described. 



In this paper we propose an algorithm for the evaluation of the hypergeometric series 
which is simpler in its practical implementation than the algorithm from [3]; the pro- 
posed algorithm is also quasi-linear in time and linear in space. The idea of working with 
a given accuracy to calculate the value of the hypergeometric series with a linear space 



complexity can be found at http:/ /numbers. computation. free.fr/Constants/Algorithms/ 
splitting.html. 

The computational complexity of the constructive real numbers and functions is 
discussed in detail in The set of constructive real numbers with quasi- linear time 
and linear space complexity of calculating their dyadic approximations will be denoted 
by Sch(FQLIN - TIME/ /LIN - SPACE) CF (CF is the Cauchy function). 

From now on, n will denote the length of the record of accuracy 2~ n of dyadic 
rational approximations. We will use \og(k) for logarithms base 2. 

1. Binary splitting method. This method is used to calculate the values of 
linearly convergent series with rational coefficients, in particular to calculate the hy- 
pergeometric series of the form 



Q _ ST- a(i) jj p(j) 

where a, b, p, and q are polynomials with integer coefficients; this series is linearly 
convergent if its partial sum 

j=0 v ; j=0 ^ KJ > 

where /i(k) is a linear function of k, differs from the exact value by not more than 2~ k : 
\S-S(ji(k))\ < 2-K 

In its classical variant the binary splitting method works as follows. Put k\ = fJ,(k). 
We consider the partial sum for some integers i\ and i 2 , < i\ < ki, < i 2 < k\, 
k < i2- 



b(i)q(ii) ■ ■ ■ q(i) 



We calculate 



P(h,h) = p(k) ■ ■ .p(i2), Q(h,i2) =q(ii)...q{ii), 

B(i 1 ,i 2 ) = b(h) ...bfa), and T(h, i 2 ) = B(i x , i 2 )Q(h, i2)S{h, 12)- 

If i\ = i 2 , then these values are calculated directly. Otherwise, the series is divided into 
two parts, left and right, and P(i\,i 2 ), Q(ii,i 2 ), and B{i\,i 2 ) are calculated for each 
part recursively. Then the values obtained are combined: 

P(h,i 2 ) = PiP r , Q (11,12) = QiQr, B(h,i 2 ) = BiB r , 
T(i 1 ,i 2 ) = B r Q r T l + B l P l T r . 

The algorithm starts with i\ = 0, i 2 = k\. After calculating T(0, ki), B(0,k\), and 
Q(0, ki), we divide T(0, k±) by B(0, k\)Q(0, k\) to get the result with the given accuracy. 
The lengths of T(0, k\) and B(0, k\)Q(Q, k\) are proportional to klog(k); therefore the 
binary splitting algorithm is quasi-linear in space; the time complexity of this algorithm 
is 0(M(k)\og(k) 2 ) 0. 



2 



2. The basic algorithm of class Sch(FQLIN - TIME/ /LIN - SPACE). 
We will modify the binary splitting method for evaluation of the hypergeometric series 
so that the algorithm is simple and is in class Sch (FQLIN — TIME / / LIN — SPACE) . 

Let's suppose that we want to calculate the values of the hypergeometric series ([1]) 
with an accuracy of 2~ n . It's enough to calculate the partial sum with accuracy 
2-(»»+l) because 

\S - S(n(n + 1))*| <\S- S(ii(n + 1))| + \S(n(n + 1)) - S(M" + !))*! 

^ 2— (n+l) _|_ 2~ ("■+!) _ 2~ n - 

here S(fi(n + 1))* is an approximate value of S(fi(n + 1)). Put r = fj,(n + 1). Take the 
minimum value k\ such that 2 kl > r; let r\ = |~^*|. We write the partial sum ([2]) as 
follows: 

P(r) = a 1 + r 2 [a 2 + r 3 [a 3 + . . . + r fcl _i[cr fcl _i + r fcl cr fcl ]]], (5) 

where 

V* o(i) tV , 

n ! 1 2 ^ a(i) A , 

50)' 2 = "^y ^ 6(0 11 ^y =6+ * 

npO') a ( 2r i) 3 ^ o(i) A p(j) , 

WY " = Wnj ^ W) 11 gy) =6+ 8 ' 



The <r£ are calculated by the usual binary splitting method for the sum ([3]), where 
i\ = (t — l)ri + 1, Z2 = i • ri — 1; for a[ we take i\ = Q, i% = r% — 1. 
We introduce the notation 

w = + 1,W = max(i, 5), 

where vl = max(|a(i)|, |6(i)|), P = max(|p(j)|, |<?(i)|) 

i=0..r j=0..r 

(here, l(u) is the length of the bit representation of u). Let's obtain an estimate of the 
computational complexity of the calculation of at and Tt- 

Lemma (1). The time complexity of the binary splitting method for calculation of at 
is bounded by 0{M(r) log(r)); the space complexity is bounded by 0(r). 

Proof. Consider an arbitrary maximal chain of recursive calls derived from calculations 
by Formulas dU). Suppose that the numbering in the chain begins with its deepest 
element: i = 1, . . . ,<;, where ? is the length of the chain, ? < [log(2ri)] . 

Using mathematical induction on i, we show that the length of the representation 
of T at level i satisfies l(TA < 2 i+1 uo + 2\ Note that /(Pi) < 2*u;, l(Qi) < 2*w, and 
/(Pj) < because there is a doubling of the length of the number representation 
when i increases. The induction begins with i = I: l(Ti) < 2u < 2 2 uj + 2 1 ; next, the 
induction step is 

l(T i+ i) < 2 1 uj + 2 1 oj + 2 i+l uj + 2 i + 1 < 2^ 1+1)+1 uj + 2 i+1 . 
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Note that, based on this inequality, T q has length 0(r) because 

Z(T f ) < C^uj < C 2 -^_l og (r) = C 2 r; 
log(r) 

here we take into account the fact that all coefficients a(i), b(i), p(j), q(j) are polyno- 
mials. 

Let's estimate the time complexity of the calculation of a' t . We must take into 
account the property of the complexity of integer multiplication: 2M(2 _1 m) < M(m) 
(quasi-linear and polynomial functions satisfy this property). Because at a tree node of 
recursive calls at level i there are C3 multiplications of 2 <; ~* +1 numbers of length at most 
2 1+1 uj + 2*, the following estimate of the number of operations required to calculate a' t 
holds: 

Time(a' t ) < C 3 J] 2^ i M(2 i+1 co + 2 l ) < C 4 ^ 2'~ i M(2 i+2 u;) 

i=l i=l 

? s 
= C 4 ^2 ? - i M(2^~( f -( i+2 »w) < C 5 ^2 ? ^2- f+ ( i+2 )M(2 ? w) 

i=l t=l 
< C 6 ?M(r) < C 7 log(r)M(r) 



(here we use the inequality for 2 ? w from the estimate of l(T^)). Final division gives 
0(M(r)) operations. 

Now we estimate the space complexity of the computation of a' t . Because at a chain 
element at level i the amount of memory used for temporary variables is Cs(2* +1 a;-|-2 J ), 
the amount of memory in all simultaneously existing recursive calls is estimated as 
follows: 

s 

Space(a' t ) < C s (2 l+l u + 2*) < C 9 (2 ? o; + 2 q ) < C w r = 0(r). 

□ 

Lemma (2). The time complexity of the binary splitting method for the calculation of 
Tt is bounded by 0(M(r) log(r)); its space complexity is bounded by 0(r). 

Proof. The estimation of the computational complexity of Tt is the same as the esti- 
mation for at due to the fact that the inequalities in the proof of lemma (1) are also 
suitable for r% (we can calculate the numerator and denominator of at using the binary 
splitting method for products). □ 

We will calculate approximate values P(r)* with accuracy 2~( n+1 ) in accordance 
with Formula ([5]) using the following iterative process: 

hi(m) = a* kl , 

hi(m) = a^_ i+1 + r^_i +2 /ii-i, i = l,...,h, (6) 
hi(m) = h,i(m) + ef, 

for i = k\ we suppose P(r)* = h^im). Here m > r (m will be chosen later), a* and 
t* are approximations of a{ and r, with accuracy 2~ m . The values hi(m) are obtained 
by discarding bits q m +iQm+2 ■ ■ ■ Qm+j of numbers hi(m) after the binary point starting 
with the (m + 1) bit: 

|£j| = \hi(m) -hi(m)\ = 0.0. . . 0q m+ iq m+2 ■ ■ ■ q m +j, (7) 
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and the sign of £j is the same as the sign of hi(m) (it is clear that [ej[ < 2 m ). 
Let's assume that the following conditions hold: 

b(i) > 2 for all i, 44 < 1 for all j. (8) 
q(j) ~ 

We prove the following two lemmas. 
Lemma (3). For every i G l..k\ 

|/ii(m)|<(t + i)nw: (9) 

Proof. We apply induction on j to hj(m). The induction base is j = 1: \hi(m)\ < 
\r\W + 2~ m < 2r{W. The induction step is (j + 1) > 2: 

\hj+i(m)\ = |cr^_ (i+1)+1 + r^ i _ (i+1)+2 /i j +£j+i\ 

< \nW + (j + i)nW + 2- m < ((j + l) + i) ri W. 

□ 

Lemma (4). The error of calculation of hk x (rn) according to scheme © is estimated 
as 

A(jfci,ro) < 2- m mkjW. 

Proof. Let's denote 

H 1 = a kl , Hi = a kl -i+i + Tk 1 -i+2H i -i, n(i,m) = \hi(m) - Hi\. 

We use the method of mathematical induction for r)(j,m) for j. The induction base is 
j = 1: 

rj(l,m) = l/fiH-ffil = |< -<7 fel | < 2~ m . 
The induction step is (j + 1) > 2: 

r?(i + l,m) = + T* kl _ {j+1)+2 hj(m) + e i+1 - 

^fci-O'+iJ+i ~ T fci-(j+i)+2#j| 

< \T*hj(m) - T v hj(m) + T v hj(m) - t v Hj \ + 2 • 2" m 

< 2- m hj(m) + T v rj(j,m) + 2 ■ 2~ m . 

Because of j9]) \hj(m)\ < (j + l)riW and, by the induction hypothesis, rj(j,m) < 
2~ m mj 2 W , we get 

r)(j + 1, m) < 2- m (j + l)nW + 2~ m mj 2 W + 2-2~ m 

< 2- m (j + l)mW + 2~ m m(j + l) 2 W = 2- m m{j + 2) 2 W. 

From A(ki,m) = n(ki,m) we get the required inequality. □ 

Lemma 4 implies that to compute P(r)* with accuracy 2~( n+1 ) it suffices to take m 
such that the following holds 

m > (n + 1) + flog(n + 1) + 2 log(r) + log(W)] , . (10) 

We denote the algorithm for the calculation of the hyper geometric series which uses 
scheme (J6j) as LinSpaceBinSplit (binary splitting method with linear space complexity). 

Algorithm LinSpaceBinSplit. The approximate value of the hypergeometric 
series. 

Input: Record of the accuracy 2~ n . 

Output: The approximate value of ([1]) with accuracy 2~ n . 
Description: 
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Table 1: Formulas for evaluation of constants 



Constant 


Formula 


Series 


e 


e = exp(lj 


exp(l) = 2 £\ =0 ^ 


TT 


tt = 16a - 4/3, 
a = arctg | , 
/3 = arctg ^9 


arctg $ 7 2 W = ^ 1 ! 1 ^ TWT ; 

arctg ( 239 ) — 2 239 Y^i=0 ( 1) 2(2i+l)239 2i 


C(3) 




V^oo (-l) I (205r i +250i+77)((i+l)!) b (i!) b 
2^i=0 2((2i+2)!) 5 



Table 2: Series for calculations of constants 



Constant 


a(i) 


b(i) 


P(j) 


q(j) 


e 


1 


2 


1 


i 


7T 


1 


2(2t + l) 


-1 


5^ 




1 


2(2* + 1) 


-1 


239 2 


C(3) 


205i 2 + 250i + 77 


2 


p(0) = 1, 


32(j + l) 5 








pO') = -i 5 





1) compute r := fj,{n + 1); choose fci so that 2 kl > r; compute r\ := [r-"|; compute 

2) find m using Formula (fTU|) : 

3) /i := cr^ (using the usual binary splitting method with accuracy 2 _m ); 

4) make loops through i from 2 to fci: 

a) calculate x>\ := o - ^ with accuracy 2~ m using the usual binary splitting 
method and V2 '■= T k 1 ~i+2 w ith accuracy 2~ m , 

b) calculate the expression h := V\ + ^2/1, 

c) assign value h rounded in accordance with (J7J) to h; 

5) write h to output. 

We estimate the time computational complexity of the algorithm taking into account 
that r and m are linearly dependent on n: 

• 0(log(n)) computations of at gives 0(M(n) log(n) 2 ); 

• 0(log(n)) computations of Tt gives 0(M(n) log(n) 2 ); 

• 0(log(n)) multiplications of numbers of the length 0{n) gives 0(M(n) log(n)); 

in total we obtain 0{M(n) log(n) 2 ) bit operations. The space complexity of algorithm 
LinSpaceBinSplit is 0(n) because in all calculations in this algorithm numbers of length 
0(n) are used. 

Theorem. The modified binary splitting algorithm for the calculation of the hypergeo- 
metric series, LinSpaceBinSplit, belongs to Sch(FQLIN - TIME/ /LIN - SPACE). 

Conclusion. In Tables [1] and [2] there are formulas and series for the calculation 
of some frequently used constants of mathematical analysis. These series converge 
linearly (see [2, 3j) and they satisfy conditions ([S]); hence, to calculate them, we can use 
LinSpaceBinSplit and therefore these constants belong to the set of constructive real 
numbers Sch(FQLLN -TIME / / LLN - SPACE) C f- The algorithm LinSpaceBinSplit 
can also be used to calculate approximations to many other constants and to the values 
of elementary functions at rational points. 
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If we use the Schonhage-Strassen algorithm for integer multiplication with a time 
complexity of 0(n log(n) log log (n)), then the time complexity of LinSpaceBinSplit will 
be 0(nlog(n) 3 loglog(n)); when we use a simple recursive method for integer multipli- 
cation with a time complexity of O(n log ^), the time complexity of LinSpaceBinSplit 
will be O(n lo s( 3 ) log(n) 2 ). 

Let's note that the series of the constant e converges with the rate 2 -0 ( nlog ( ra )) ) so the 
time complexity of the calculation of e using LinSpaceBinSplit will be 0(M{n) log(n)) 

and its space complexity will be O ( , 5 , ) . 
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